/******************************************************************************
* File: main_plaquette.cpp
*
* Author: Vishal Rao
* Created: 02/08/24
* Description: This program contains the function to calculate the wilson loop with hyp smearing and generation of more points using bresenham algorithm described in :https://arxiv.org/pdf/hep-lat/0005018.pdf
*****************************************************************************
* main_plaquette.cpp
* This is just an example how a very very basic program works. Look at src/testing/main_GeneralOperatorTest.cpp to see
* how to write more advanced GPU code.
*
*/
#include "../simulateqcd.h"
/* #include "../modules/hyp/hypParameters.h" was already included in ../modules/hypSmearing.h */
/* #include "../modules/hyp/smearParameters.h" */
#include "../modules/hyp/hypSmearing.h"
/* #include "../modules/hyp/hypSmearing.cpp" */
#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#define PREC double
template<class floatT,size_t HaloDepth>
struct CalcWilson{
//Gauge accessor to access the gauge field
SU3Accessor<floatT> SU3Accessor;
int _length=1;
int _height=1;
int _nu_in=0; // nu_in is the direction of fictitious time , in this direction we don't do smearing.
//Constructor to initialize all necessary members.
CalcWilson(Gaugefield<floatT,true,HaloDepth> &gauge, int &length, int &height, int &nu_in) : SU3Accessor(gauge.getAccessor()), _length(length),_height(height), _nu_in(nu_in){
}
__device__ __host__ SU3<floatT> HalfLoop(int &length, int &height,gSite site,int &mu, int &nu)
{
// nu_in is the direction of fictitious time , in this direction we don't do smearing. // we have not taken site variable as reference as if it were take as reference, it will update the values at same mEmory location
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1); //initialize to identity matrix
typedef GIndexer<All, HaloDepth> GInd;
for (int il=0;il<length;il++){
temp *= SU3Accessor.getLink(GInd::getSiteMu(site, mu));
site=GInd::site_up(site, mu);
}
for (int ih=0;ih<height;ih++){
temp *= SU3Accessor.getLink(GInd::getSiteMu(site, nu));
site=GInd::site_up(site,nu);
}
return temp;
}
//This is the operator that is called inside the Kernel // nu_in is the direction of fictitious time , in this direction we don't do smearing.
__device__ __host__ floatT operator()(gSite site) {
/// We need to choose the type of indexer. The first template is the layout of the lattice. // nu_in is the direction of fictitious time , in this direction we don't do smearing.
typedef GIndexer<All, HaloDepth> GInd;
int l=_length;
int h=_height;
/// Define a SU(3) matrix
SU3<floatT> temp1,temp2;
floatT result = 0;
int mu1,nu1,mu2,nu2;
switch (_nu_in) {
case 0: {nu1=0;nu2=0;mu1=1;mu2=2;}
break;
case 1: {nu1=1;nu2=1;mu1=0;mu2=2;}
break;
case 2: {nu1=2;nu2=2;mu1=0;mu2=1;}
break;
}
temp1=HalfLoop(l,h,site,mu1,nu1);
temp2=HalfLoop(h,l,site,nu1,mu1);
result+=tr_d(temp2*dagger(temp1));
temp1=HalfLoop(l,h,site,mu2,nu2);
temp2=HalfLoop(h,l,site,nu2,mu2);
result+=tr_d(temp2*dagger(temp1));
return result/2.0;
}
};
template<class floatT,size_t HaloDepth>
struct CalcWilsonBresenham{
int _x,_y,_z, _ex_dir;
SU3Accessor<floatT> SU3Accessor;
CalcWilsonBresenham(Gaugefield<floatT,true,HaloDepth> &gauge, int &x,int &y, int &z, int &ex_dir) : SU3Accessor(gauge.getAccessor()), _x(x),_y(y),_z(z), _ex_dir(ex_dir){}
//This is the operator that is called inside the Kernel
__device__ __host__ SU3<floatT> BresenhamLoopXYL(int &x,int &y, int &z, gSite site, int &mu, int &nu)
{
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1); //initialize to identity matrix
int cmax,cmin,cmax_dir,cmin_dir;
if (x>=y) {cmax=x; cmin=y; cmax_dir=mu; cmin_dir=nu;}
else {cmax=y; cmin=x; cmax_dir=nu; cmin_dir=mu;}
typedef GIndexer<All, HaloDepth> GInd;
int cmax2=2*cmax;
int cmin2=2*cmin;
int chi=cmin2-cmax;
for (int ix=0;ix<cmax;ix++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmax_dir));//step in cmax
site=GInd::site_up(site,cmax_dir);
if(chi>=0)
{
chi-=cmax2;
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmin_dir)); //step in cmin
site=GInd::site_up(site,cmin_dir);
}
chi+=cmin2;
}
for(int iz=0;iz<z;iz++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,3-mu-nu)); //step in cmin
site=GInd::site_up(site,3-mu-nu);
}
return temp;
}
__device__ __host__ SU3<floatT> BresenhamLoopXYU(int &x,int &y, int &z, gSite site, int &mu, int &nu)
{
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1); //initialize to identity matrix
int cmax,cmin,cmax_dir,cmin_dir;
if (x>=y) {cmax=x; cmin=y; cmax_dir=mu; cmin_dir=nu;}
else {cmax=y; cmin=x; cmax_dir=nu; cmin_dir=mu;}
typedef GIndexer<All, HaloDepth> GInd;
int cmax2=2*cmax;
int cmin2=2*cmin;
int chi=cmin2-cmax;
for(int iz=0;iz<z;iz++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,3-mu-nu)); //step in cmin
site=GInd::site_up(site,3-mu-nu);
}
for (int ix=0;ix<cmax;ix++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmax_dir));//step in cmax
site=GInd::site_up(site,cmax_dir);
if(chi>=0)
{
chi-=cmax2;
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmin_dir)); //step in cmin
site=GInd::site_up(site,cmin_dir);
}
chi+=cmin2;
}
return temp;
}
__device__ __host__ floatT operator()(gSite site) {
typedef GIndexer<All, HaloDepth> GInd;
/// Define a SU(3) matrix
floatT result = 0;
SU3<floatT> temp1,temp2;
int mu,nu;
switch (_ex_dir) {
case 0: mu=1;nu=2;break;
case 1: mu=0;nu=2;break;
case 2: mu=0;nu=1;break;
}
temp1=BresenhamLoopXYL(_x,_y,_z,site,mu,nu);
temp2=BresenhamLoopXYU(_x,_y,_z,site,mu,nu);
result+=tr_d(temp1*dagger(temp2));
//replacing x<->y
temp1=BresenhamLoopXYL(_y,_x,_z,site,mu,nu);
temp2=BresenhamLoopXYU(_y,_x,_z,site,mu,nu);
result+=tr_d(temp1*dagger(temp2));
return result/2.0;
}
};
//Function to compute the wilson loop using the above struct CalcWilson.
template<class floatT, size_t HaloDepth>
floatT WilsonLoop(Gaugefield<floatT,true, HaloDepth> &gauge, LatticeContainer<true,floatT> &redBase, int &l, int &h, int &ex_dir){
typedef GIndexer<All,HaloDepth> GInd;
const size_t elems = GInd::getLatData().vol4;
//Make sure, redBase is large enough
redBase.adjustSize(elems);
// we will iterate the process of finding wilson loop over whole bulk.
redBase.template iterateOverBulk<All, HaloDepth>(CalcWilson<floatT, HaloDepth>(gauge, l,h,ex_dir)); //ex_dir is excluded direction for smearing
//Do the final reduction
floatT Wloop;
redBase.reduce(Wloop, elems);
//Normalize the result
Wloop /= (GInd::getLatData().globalLattice().mult()*3); //3 colors
return Wloop;
}
//defining function for bresenham
template<class floatT, size_t HaloDepth>
floatT WilsonLoopBresenham(Gaugefield<floatT,true, HaloDepth> &gauge, LatticeContainer<true,floatT> &redBase, int &x, int &y,int &z, int &ex_dir){
typedef GIndexer<All,HaloDepth> GInd;
const size_t elems = GInd::getLatData().vol4;
//Make sure, redBase is large enough
redBase.adjustSize(elems);
// we will iterate the process of finding wilson loop over whole bulk.
redBase.template iterateOverBulk<All, HaloDepth>(CalcWilsonBresenham<floatT, HaloDepth>(gauge,x,y,z, ex_dir));
//Do the final reduction
floatT Wloop;
redBase.reduce(Wloop, elems);
//Normalize the result
const int n_colors=3;
Wloop /= (GInd::getLatData().globalLattice().mult()*n_colors); //3 loops(one averaged loop per plane) times 3 colors=9.
return Wloop;
}
//
// I am going to paste the hypSmearing.cpp file content here for the moment.
//HYP smearing
// also call updateAll()
template<class floatT, bool onDevice, size_t HaloDepth, CompressionType comp>
void HypSmearing<floatT, onDevice, HaloDepth, comp>::Su3Unitarize(Gaugefield<floatT, onDevice, HaloDepth, comp> &gauge_out, Gaugefield<floatT, onDevice, HaloDepth, comp> &gauge_base){
//4 is link number
HypStaple<floatT, HaloDepth, comp, 4> su_3_unitarize(gauge_out.getAccessor(), gauge_base.getAccessor(), gauge_base.getAccessor(), gauge_base.getAccessor());
gauge_out.iterateOverBulkAllMu(su_3_unitarize);
if(update_all)gauge_out.updateAll();
}
template<class floatT, bool onDevice, size_t HaloDepth, CompressionType comp>
void HypSmearing<floatT, onDevice, HaloDepth, comp>::SmearAll(Gaugefield<floatT, onDevice, HaloDepth, comp> &gauge_out)
{
// create level 1 fields
_dummy.iterateOverBulkAllMu(staple3_lvl1_10);
_gauge_lvl1_10 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_10, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_20);
_gauge_lvl1_20 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_20, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_30);
_gauge_lvl1_30 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_30, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_21);
_gauge_lvl1_21 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_21, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_31);
_gauge_lvl1_31 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_31, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_32);
_gauge_lvl1_32 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_32, _gauge_base);
// now that we have level 1 fields, create level 2 staples
// note: the order of the gauge fields goes in ascending order (10 < 20 < 30, 10 < 21 < 31, 20 < 21 < 32, 30 < 31 < 32)
// this is ASSUMED by HypStaple<floatT, HaloDepth, comp, 2>; DO NOT change this order without also modifying HypStaple<floatT, HaloDepth, comp, 2> and threeLinkStaple_second_level
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_0(_gauge_lvl1_10.getAccessor(), _gauge_lvl1_20.getAccessor(), _gauge_lvl1_30.getAccessor(), _dummy.getAccessor(), 0);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_1(_gauge_lvl1_10.getAccessor(), _gauge_lvl1_21.getAccessor(), _gauge_lvl1_31.getAccessor(), _dummy.getAccessor(), 1);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_2(_gauge_lvl1_20.getAccessor(), _gauge_lvl1_21.getAccessor(), _gauge_lvl1_32.getAccessor(), _dummy.getAccessor(), 2);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_3(_gauge_lvl1_30.getAccessor(), _gauge_lvl1_31.getAccessor(), _gauge_lvl1_32.getAccessor(), _dummy.getAccessor(), 3);
//second level fields
_dummy.iterateOverBulkAllMu(staple3_lvl2_0);
_gauge_lvl2_0 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_0, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl2_1);
_gauge_lvl2_1 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_1, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl2_2);
_gauge_lvl2_2 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_2, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl2_3);
_gauge_lvl2_3 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_3, _gauge_base);
// now that we have level 2 fields, create level 3 staple
HypStaple<floatT, HaloDepth, comp, 1> staple3_lvl3(_gauge_lvl2_0.getAccessor(), _gauge_lvl2_1.getAccessor(), _gauge_lvl2_2.getAccessor(), _gauge_lvl2_3.getAccessor());
_dummy.iterateOverBulkAllMu(staple3_lvl3);
// OLD VERSION, (MAYBE) DOES NOT WORK FOR SOME REASON
//_gauge_lvl2_0 = (1-params.alpha_1) * _gauge_base + params.alpha_1/6 * _dummy; //reused _gauge_lvl2_0
//Su3Unitarize(_gauge_lvl2_0);
// NEW VERSION, USES EXTRA FIELD RATHER THAN REUSE _gauge_lvl2_0
gauge_out = (1-params.alpha_1) * _gauge_base + params.alpha_1/6 * _dummy; //reused _gauge_lvl2_0
Su3Unitarize(gauge_out, _gauge_base);
}
template<class floatT, bool onDevice, size_t HaloDepth, CompressionType comp>
void HypSmearing<floatT, onDevice, HaloDepth, comp>::SmearAll_spatial(Gaugefield<floatT, onDevice, HaloDepth, comp> &gauge_out, int &ex_dir)
{
switch (ex_dir) {
case 0:
{ _dummy.iterateOverBulkAllMu(staple3_lvl1_10);
_gauge_lvl1_10 = (1-params.alpha_2) * _gauge_base + params.alpha_2/2 * _dummy;
Su3Unitarize(_gauge_lvl1_10, _gauge_base); //_gauge_level1_10 is su3 now
_dummy.iterateOverBulkAllMu(staple3_lvl1_20);
_gauge_lvl1_20 = (1-params.alpha_2) * _gauge_base + params.alpha_2/2 * _dummy;
Su3Unitarize(_gauge_lvl1_20, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_30);
_gauge_lvl1_30 = (1-params.alpha_2) * _gauge_base + params.alpha_2/2 * _dummy;
Su3Unitarize(_gauge_lvl1_30, _gauge_base);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_0(_gauge_lvl1_10.getAccessor(), _gauge_lvl1_20.getAccessor(), _gauge_lvl1_30.getAccessor(), _dummy.getAccessor(), 0);
_dummy.iterateOverBulkAllMu(staple3_lvl2_0);
gauge_out = (1-params.alpha_1) * _gauge_base + params.alpha_1/4 * _dummy; //reused _gauge_lvl2_0
Su3Unitarize(gauge_out, _gauge_base);
} break;
case 1:
{ _dummy.iterateOverBulkAllMu(staple3_lvl1_10);
_gauge_lvl1_10 = (1-params.alpha_2) * _gauge_base + params.alpha_2/2 * _dummy;
Su3Unitarize(_gauge_lvl1_10, _gauge_base); //_gauge_level1_10 is su3 now
_dummy.iterateOverBulkAllMu(staple3_lvl1_21);
_gauge_lvl1_21 = (1-params.alpha_2) * _gauge_base + params.alpha_2/2 * _dummy;
Su3Unitarize(_gauge_lvl1_21, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_31);
_gauge_lvl1_31 = (1-params.alpha_2) * _gauge_base + params.alpha_2/2 * _dummy;
Su3Unitarize(_gauge_lvl1_31, _gauge_base);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_1(_gauge_lvl1_10.getAccessor(), _gauge_lvl1_21.getAccessor(), _gauge_lvl1_31.getAccessor(), _dummy.getAccessor(), 1);
_dummy.iterateOverBulkAllMu(staple3_lvl2_1);
gauge_out = (1-params.alpha_1) * _gauge_base + params.alpha_1/4 * _dummy; //reused _gauge_lvl2_0
Su3Unitarize(gauge_out, _gauge_base);
} break;
case 2:
{ _dummy.iterateOverBulkAllMu(staple3_lvl1_20);
_gauge_lvl1_20 = (1-params.alpha_2) * _gauge_base + params.alpha_2/2 * _dummy;
Su3Unitarize(_gauge_lvl1_20, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_21);
_gauge_lvl1_21 = (1-params.alpha_2) * _gauge_base + params.alpha_2/2 * _dummy;
Su3Unitarize(_gauge_lvl1_21, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_32);
_gauge_lvl1_32 = (1-params.alpha_2) * _gauge_base + params.alpha_2/2 * _dummy;
Su3Unitarize(_gauge_lvl1_32, _gauge_base);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_2(_gauge_lvl1_20.getAccessor(), _gauge_lvl1_21.getAccessor(), _gauge_lvl1_32.getAccessor(), _dummy.getAccessor(), 2);
_dummy.iterateOverBulkAllMu(staple3_lvl2_2);
gauge_out = (1-params.alpha_1) * _gauge_base + params.alpha_1/4 * _dummy; //reused _gauge_lvl2_0
Su3Unitarize(gauge_out, _gauge_base);
} break;
}
}
#define CLASS_INIT(floatT,HALO) \
template class HypSmearing<floatT,true,HALO,R18>;
INIT_PH(CLASS_INIT)
int main(int argc, char *argv[]) {
stdLogger.setVerbosity(DEBUG);
/// Initialize parameter class. This class can also read parameter from textfiles!
LatticeParameters param;
/// Initialize the Lattice dimension
const int LatDim[] = {32, 32, 32, 8};
const int NodeDim[] = {1, 1, 1, 1};
/// Just pass these dimensions to the parameter class
param.latDim.set(LatDim);
param.nodeDim.set(NodeDim);
/// Initialize a timer
StopWatch<true> timer;
/// Initialize the CommunicationBase. This class handles the communitation between different Cores/GPU's.
CommunicationBase commBase(&argc, &argv, true);
commBase.init(param.nodeDim());
const size_t HaloDepth = 1; //since it is not multi-gpu code
/// highlight the output differently.
rootLogger.info("Initialize Lattice");
/// Initialize the Indexer on GPU and CPU.
initIndexer(HaloDepth,param,commBase);
typedef GIndexer<All,HaloDepth> GInd;
rootLogger.info("Initialize Gaugefield");
Gaugefield<PREC, true,HaloDepth> gauge(commBase);
typedef double floatT;
GaugeAction<floatT, true, HaloDepth, R18> gaugeaction(gauge);
/// Initialize gaugefield with unity-matrices.
gauge.one();
/// Initialize LatticeContainer. This is in principle the "array", where the values of the plaquette are
/// stored which are summed up in the end
LatticeContainer<true,PREC> redBase(commBase);
/// We need to tell the Reductionbase how large our Array will be
redBase.adjustSize(GInd::getLatData().vol4);
int eqm_point=155;
int last_point=1150;
int node;
sscanf(argv[1], "%d", &node);
//sscanf(argv[2], "%d", &height);
int no_smearing_steps=1;
for (int i=eqm_point; i<=last_point;i+=5)
{
rootLogger.info("Read configuration");
gauge.readconf_nersc("/root/project1/build_SIMULATeQCD/applications/try_20_output_dir/after_eqm/node"+std::to_string(node)+"/l328f21b6285m0039185m0783706a_"+std::to_string(node)+"."+std::to_string(i));
gauge.updateAll();
Gaugefield<PREC, true,HaloDepth> gauge_out(commBase);
gauge_out.one(); //initialize the gauge field variable
Gaugefield<PREC, true,HaloDepth> gauge_copy(commBase);
gauge_copy.one(); //initialize the gauge field variable
gauge_copy=gauge;
for(int ex_dir=0;ex_dir<3;ex_dir++)
{
for(int j=1;j<=no_smearing_steps;j++)
{
HypSmearing<floatT, true, HaloDepth, R18> hypsmearing(gauge); //smearing of gauge
/* hypsmearing.SmearAll(gauge_out); //and storing at gauge_out */
hypsmearing.SmearAll_spatial(gauge_out, ex_dir); //and storing at gauge_out
gauge=gauge_out;
}
gauge=gauge_copy;
const int Ns= LatDim[0];
for (int length=1; length<=1; length++)
{
for (int height=1; height<Ns; height++)
{
PREC Wloop = 0;
/* timer.start(); */
/// compute wilson loop with smeared gauge_out
Wloop = WilsonLoop<PREC,HaloDepth>(gauge_out, redBase, length,height, ex_dir );
printf("Node:%d,Smear_count:%d,Length=%d,Height=%d,ex_dir=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,length, height,ex_dir,Wloop);
/* rootLogger.info("Reduced Rectangle from rhmc : " , gaugeaction.rectangle()); */
/* rootLogger.info("Reduced plaquette from rhmc : " , gaugeaction.plaquette()); */
}
}
// we shall now implement bresenham algorithm to calculate the new wilson loops.int x_coord=length;
/* for (int a=1;a<Ns;a++) */
/* { */
/* for(int b=a;b<Ns;b++) */
/* { */
/* for(int c=1;c<Ns;c++) */
/* { */
/* PREC Wloop=WilsonLoopBresenham<PREC,HaloDepth>(gauge_out, redBase,a,b,c, ex_dir ); */
/* printf("Node:%d,Smear_count:%d,Length=%f,Height=%d,ex_dir=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,pow((a*a+b*b),0.5), c,ex_dir,Wloop); */
/* } */
/* } */
/* } */
}
}
return 0;
}